#include "allincludes.h"
#include <vector>
using namespace std;
int OneLine::NormalizeData()
{
	if(!bIsNormalized)
	{
		/*int USEABLESTARTINDEX;
		if(binsLengthX==1)
			USEABLESTARTINDEX=USEABLESTARTINDEX_7_5M;
		else if(binsLengthX==2)
			USEABLESTARTINDEX=USEABLESTARTINDEX_15M;
		else
		{
			cout<<"Can't handle this resolution!"<<endl;
			exit(FAILED);
		}*/
		for(int i=USEABLESTARTINDEX;i<binsNum;i++)
		{
			data[i] = data[i]-analogData[BACKGROUNDINDEX]+ini.farrestAfterpulseAverage*analogData[ENERGYINDEX];
			if(i<=afterpulseEndIndex)
			{
				data[i]=data[i]-afterpulse[i]*analogData[ENERGYINDEX];
			}
			if(i<=overlapEndIndex)
			{
				data[i]=data[i] * overlapCoefficient[i];
			}
			data[i]=data[i]*pow(IndexToDistance(i),2);

			//20080318 From now on, the Normalized data is devided by LidarConst and Energy.
			data[i]=data[i]/analogData[ENERGYINDEX]/ini.lidarConst;
			//20080322 Set the minus data to zero...
			if(data[i]<0.0)
				data[i]=0.0;

		}
		if(USEABLESTARTINDEX>FIRSTDATAINDEX)          //Make the first several un-useable data the same as the first useable data.
		{
			double aver = 0.0;
			for(int i=USEABLESTARTINDEX; i < USEABLESTARTINDEX + 5; ++i)
			{
				aver += data[i];
			}
			aver /= 5.0;

			for(int i=FIRSTDATAINDEX;i<USEABLESTARTINDEX;i++)
			{
				data[i]=aver;
			}
		}
		bIsNormalized=true;
	}
	return SUCCEED;
}
OneLine OneLine::GetSigmaA(int neededHeightIndex, float newSa,float newLidarConst )
{
	if(bIsSigma)
	{
		return *this;
	}
	else
	{
		if(!bIsNormalized)
		{
			NormalizeData();
		}
	/*	int FIRSTDATAINDEX;
		if(binsLengthX==1) FIRSTDATAINDEX=FIRSTDATAINDEX_7_5M;
		else if(binsLengthX==2) FIRSTDATAINDEX=FIRSTDATAINDEX_15M;
		else {cout<<"Can't handle this resolution!"<<endl;exit(FAILED);}*/
		/*int USEABLESTARTINDEX;
		if(binsLengthX==1)
			USEABLESTARTINDEX=USEABLESTARTINDEX_7_5M;
		else if(binsLengthX==2)
			USEABLESTARTINDEX=USEABLESTARTINDEX_15M;
		else
		{
			cout<<"Can't handle this resolution!"<<endl;
			exit(FAILED);
		}*/

		if(neededHeightIndex<=0)
			neededHeightIndex=binsNum-1;
		if(newLidarConst<=0.1)
			newLidarConst=ini.lidarConst;
		if(newSa<=0.1)
			newSa=ini.sa;
		OneLine sigmaALine(*this);
		
		vector<float> exp_term(neededHeightIndex+1);
		vector<float> old_x_exp(neededHeightIndex+1);
		vector<float> cal(neededHeightIndex+1);

		for(int i=FIRSTDATAINDEX;i<=neededHeightIndex && i<binsNum;i++)
		{
			exp_term[i] = exp( -2.0 * (newSa - SM) * integratedBetaM[i]);
			old_x_exp[i] = data[i] * exp_term[i];
			if(i == FIRSTDATAINDEX)
				cal[i] = old_x_exp[i] * 0.0075 * binsLengthX * 0.5;
			else
				cal[i] = cal[i - 1] + old_x_exp[i - 1] * 0.0075 * binsLengthX * 0.5 + old_x_exp[i] * 0.0075 * binsLengthX * 0.5;

//			for(int j=FIRSTDATAINDEX;j<=i;j++)  //Maybe we should manage to fill up the gap between first data and first usable data
//			{
//				cal += data[j]*exp(-2*(newSa-SM)*integratedBetaM[j]);
//			}
//			cal *= 0.0075*binsLengthX;
			//20080318 From now on, the Normalized data is devided by LidarConst and Energy.
			sigmaALine.data[i]=( old_x_exp[i] / ( 1.0 - 2*newSa*cal[i] ) - betaM[i] )*newSa;
			//20080322 We want no minux data...
			if(sigmaALine.data[i]<0.0)
				sigmaALine.data[i]=0.0;
		}

		sigmaALine.bIsSigma=true;
		sigmaALine.maxIndexOfSigmaedData=neededHeightIndex;
		return sigmaALine;
	}
}

/*int OneLine::DoRoutine()
{
	NormalizeData();
	///
	///
	
	return SUCCEED;
}*/



float OneLine::GetSa(float AOD, int neededHeightIndex, string outFileName, float minSa , float maxSa )
{
	if(bIsSigma)
	{
		cout<<"Can't get Sa from Sigma data... , using old Sa"<<endl;
		return ini.sa;
	}
	else
	{
		if(bIsNormalized)
		{
			NormalizeData();
		}
		if (outFileName=="")
		{
			outFileName=ini.outputPath+ GetTimeString() +"_"+ GetAverageMinString()+"min.average.Sa_AOD.txt";
		}
		ofstream outFile(outFileName.c_str());

		if(!outFile)
		{
			cout<<"Can't create or open file : "<<outFileName<<endl;
			exit(FAILED);
		}
		
		float closestAOD=-9999;
		float resultSa=ini.sa;

		outFile<<setw(16)<<"new Sa"<<setw(16)<<"new AOD"<<"  --and Target AOD is "<<AOD<<endl;
		cout<<setw(12)<<"Result Sa"<<setw(12)<<"new Sa"<<setw(12)<<"new AOD"<<setw(12)<<"Target AOD"<<endl;
		for(float newSa=minSa; newSa<=maxSa;newSa+=0.1)
		{
			OneLine temp=GetSigmaA(neededHeightIndex, newSa);
			float newAOD=temp.GetAOD(neededHeightIndex, newSa);

			outFile<<setw(16)<<setprecision(6)<<newSa<<setw(16)<<setprecision(6)<<newAOD<<endl;

			if( fabs(newAOD-AOD) <= fabs(closestAOD-AOD) )
			{
				closestAOD=newAOD;
				resultSa=newSa;
			}
			cout<<setw(12)<<setprecision(3)<<resultSa<<setw(12)<<setprecision(3)<<newSa<<setw(12)<<newAOD<<setw(12)<<AOD<<endl;

			if(newAOD>10 ||newAOD<-10)
			{
				cout<<"  -- got panic... breaking down..."<<endl<<endl;
				break;
			}
		}
		if ( fabs(closestAOD-AOD)>0.2 )
		{
			cout<<"Can't find proper sa, using previous one"<<endl;
			outFile<<"Can't find proper sa, using previous one : "<<ini.sa<<endl;
			outFile.close();

			return ini.sa;
		}
		else
		{
			cout<<"Result sa: "<<resultSa<<"   observed AOD: "<<AOD<<"   calculated AOD: "<<closestAOD<<endl;
			outFile<<"Result sa: "<<resultSa<<"   observed AOD: "<<AOD<<"   calculated AOD: "<<closestAOD<<endl;
			outFile.close();
			
			return resultSa;
		}
	}

}

float OneLine::GetAOD(int neededHeightIndex, float newSa,float newLidarConst)
{
	/*int USEABLESTARTINDEX;
	if(binsLengthX==1)
		USEABLESTARTINDEX=USEABLESTARTINDEX_7_5M;
	else if(binsLengthX==2)
		USEABLESTARTINDEX=USEABLESTARTINDEX_15M;
	else
	{
		cout<<"Can't handle this resolution!"<<endl;
		exit(FAILED);
	}*/

	float nowSa, nowLidarConst;
	nowSa=newSa;
	nowLidarConst=newLidarConst;
	if(neededHeightIndex<=0)
		neededHeightIndex=binsNum;
	if(newSa<=0.1)
		nowSa=ini.sa;
	if(newLidarConst<=0.1)
		nowLidarConst=ini.lidarConst;

	OneLine tempLine;
	if(bIsSigma)
	{
		if(neededHeightIndex>maxIndexOfSigmaedData)
			cout<<"  --Warning: you might be using GetAOD() with an incompleted sigma_a line..."<<endl;
		tempLine=*this;
	}
	else
	{
		tempLine=GetSigmaA(neededHeightIndex, nowSa, nowLidarConst);
	}
	float nowAOD=0.0;
	for(int i=FIRSTDATAINDEX;i<=neededHeightIndex && i< binsNum;i++)
	{
		nowAOD += tempLine.data[i];
	}
	nowAOD *= 0.0075*binsLengthX*SINALPHA; //??

	return nowAOD;

	
}
int OneLine::GetOverlapCoefficient(float linearStartDistance, float linearEndDistance,string outFileName)
{
	string outFileName2=outFileName+".Detail";
	if(bIsSigma||bIsNormalized)
	{
		cout<<"Can't get Overlap Coefficient from Sigma data or Normalized data... "<<endl;
		return FAILED;
	}
	else
	{
		OneLine temp(*this);
		//OneLine temp2;

		for(int i=USEABLESTARTINDEX;i<temp.binsNum;i++)
		{
			temp.data[i] = temp.data[i]-temp.analogData[BACKGROUNDINDEX]+temp.ini.farrestAfterpulseAverage*temp.analogData[ENERGYINDEX];
			if(i<=afterpulseEndIndex)
			{
				temp.data[i]=temp.data[i]-afterpulse[i]*temp.analogData[ENERGYINDEX];
			}
			temp.data[i]=temp.data[i]*pow(IndexToDistance(i),2);
		//	temp2=temp;

			if (temp.data[i]<LOGLOWLIMIT)
				temp.data[i]=LOGLOWLIMIT;      
			temp.data[i]=log(temp.data[i]);
		}

		int nearIndex = DistanceToIndex(linearStartDistance);
		int farIndex =DistanceToIndex(linearEndDistance);
		if(farIndex>=temp.binsNum)
			farIndex=temp.binsNum-1;
		if(nearIndex<10)			//
			nearIndex=10;   //
		if(farIndex<nearIndex)
			farIndex=nearIndex;

		float tempFar=IndexToDistance(farIndex);
		float tempNear=IndexToDistance(nearIndex);

		//linear fit.... -__- |||b
		float xBar=0;
		float yBar=0;
		float lxx=0;
		float lxy=0;
		float lyy=0;
		float rxy=0;
		for(int i=nearIndex; i<=farIndex;i++)
		{
			xBar += IndexToDistance(i);
			yBar += temp.data[i];
		}
		xBar /= (farIndex-nearIndex+1);
		yBar /= (farIndex-nearIndex+1);

		for(int i=nearIndex; i<=farIndex;i++)
		{
			lxx += pow((IndexToDistance(i)-xBar),2);
			lxy +=  (IndexToDistance(i)-xBar)*(temp.data[i]-yBar);
			lyy += pow( (temp.data[i]-yBar), 2);
		}
		
		rxy=lxy/sqrt(lxx*lyy);

		cout<<"rxy= "<<rxy<<endl;
		if(fabs(rxy)<0.9)
		{
			cout<<"  -- which means you can hardly trust this linear fit result."<<endl;
		}
		float b=lxy/lxx;
		float a=yBar-b*xBar;
		cout<<"y=a+bx, a= "<<a<<"   , b= "<<b<<endl;
		cout<<"Visible Distance is about: "<< 3.91*2/fabs(b)<<" km."<<endl;
		//float *list;
		//list= new float[highIndex-lowIndex+1];
		//
		//float sum=0.0;
		//for(int i=lowIndex;i<=highIndex;i++)
		//{
		//	list[i-lowIndex]=(float)(data[i]/(analogData[ENERGYINDEX]*betaM[i]*exp(-2*AOD-2*integratedBetaM[i]*SM)));
		//	sum += list[i-lowIndex];
		//}
		//float result=sum/(highIndex-lowIndex+1);
		//cout<<"The Lidar Const is : "<<result<<endl;
		if (outFileName=="")
		{
			outFileName=ini.outputPath + GetTimeString() +"_"+ GetAverageMinString()+"min.OverlapCoef.txt";
			outFileName2=ini.outputPath + GetTimeString() +"_"+ GetAverageMinString()+"min.OverlapCoef.Detail.txt";
		}
		ofstream outFile(outFileName.c_str());
		ofstream outFile2(outFileName2.c_str());

		if(!outFile)
		{
			cout<<"Can't create or open file : "<<outFileName<<endl;
			return FAILED;
		}
		else if(!outFile2)
		{
			cout<<"Can't create or open file : "<<outFileName2<<endl;
			return FAILED;
		}
		else
		{		
			outFile<<"This is an Overlap Coef data file , made with Lidar record"<<temp.GetTimeString()<<" , with "<<temp.GetAverageMinString()<<" mins average"<<endl;
			outFile<<"The Linear part of the log(data) curve starts from "<<tempNear<<" km, ends at "<<tempFar<<" km, with the equation : y= "<<a<<" [+] "<<b<<" x, rxy: "<<rxy<<endl;
			outFile<<"Please find out where the Overlap Coef is 1, and trunc this file to make it suitable for the program's input. and edit the following data count number... p.s. , there should be 3 comment lines left..."<<endl;
			outFile<<temp.binsNum-USEABLESTARTINDEX<<"           "<<temp.binsLengthX<<endl;
			
			outFile2<<"This is an Overlap Coef Detail file , made with Lidar record"<<temp.GetTimeString()<<" , with "<<temp.GetAverageMinString()<<" mins average"<<endl;
			outFile2<<"The Linear part of the log(data) curve starts from "<<tempNear<<" km, ends at "<<tempFar<<" km, with the equation : y= "<<a<<" [+] "<<b<<" x, rxy: "<<rxy<<endl;
			outFile2<<"This Data is for further process. Don't use it as Overlap Coefficient directly"<<endl;
			outFile2<<temp.binsNum-USEABLESTARTINDEX<<"           "<<temp.binsLengthX<<endl;
			
			for(int i=USEABLESTARTINDEX;i<temp.binsNum;i++)
			{
				outFile<<setw(8)<<setprecision(6)<<i<<setw(16)<<setprecision(6)<<IndexToDistance(i)<<setw(16)<<setprecision(6)<<1.0/exp(temp.data[i]-a-b*IndexToDistance(i))<<endl;
				outFile2<<setw(8)<<setprecision(6)<<i<<setw(16)<<setprecision(6)<<IndexToDistance(i)<<setw(16)<<setprecision(6)<<temp.data[i]<<endl;
				
				/*if(fabs(temp2.data[i])<0.00001)
					outFile<<setw(8)<<setprecision(6)<<i<<setw(16)<<setprecision(6)<<IndexToDistance(i)<<setw(16)<<setprecision(6)<<0.0<<endl;
				else
					outFile<<setw(8)<<setprecision(6)<<i<<setw(16)<<setprecision(6)<<IndexToDistance(i)<<setw(16)<<setprecision(6)<<exp(a+b*IndexToDistance(i))/temp2.data[i]<<endl;*/
			}
			outFile.close();
			outFile2.close();
			return SUCCEED;
		}
	}
}


float OneLine::GetLidarConst(float AOD, float lowHeight, float highHeight, bool bOutputDetail,string outFileName)
{
	if(bIsSigma)
	{
		cout<<"Can't get Lidar Const from Sigma data... , returning old LC"<<endl;
		return ini.lidarConst;
	}
	else
	{
		if(!bIsNormalized)
		{
			NormalizeData();
		}

		int lowIndex = HeightToIndex(lowHeight);
		int highIndex =HeightToIndex(highHeight);
		
		if(highIndex>=binsNum)
			highIndex=binsNum-1;
		if(lowIndex<10)			//
			lowIndex=10;   //
		if(highIndex<lowIndex)
			highIndex=lowIndex;

		float tempHigh=IndexToHeight(highIndex);
		float tempLow=IndexToHeight(lowIndex);

		float *list;
		list= new float[highIndex-lowIndex+1];
		
		cout<<setw(12)<<"LidarConst"<<"   "<<setw(8)<<"index"<<"   "<<setw(12)<<"data"<<endl;
		float sum=0.0;
		for(int i=lowIndex;i<=highIndex;i++)
		{
			//080322 removing the energy devided here because of change of the defination of normalising;
			list[i-lowIndex]=(float)(data[i]/(betaM[i]*exp(-2*AOD-2*integratedBetaM[i]*SM))*ini.lidarConst);
			if(data[i]<0)
				cout<<"! ";
			else
				cout<<"  ";
			cout<<setw(10)<<setprecision(5)<<list[i-lowIndex]<<"   "<<setw(8)<<i<<"   "<<setw(12)<<data[i]<<endl;

			sum += list[i-lowIndex];
		}
		float result=sum/(highIndex-lowIndex+1);
		cout<<"The Lidar Const Calculated and Averaged is : "<<result<<endl;



		if (outFileName=="")
		{
			outFileName=ini.outputPath+  GetTimeString() +"_"+ GetAverageMinString()+"min.average.LidarConst.txt";
		}
		ofstream outFile(outFileName.c_str(),ios::app);

		if(!outFile)
		{
			cout<<"Error : Can't create or open file : "<<outFileName<<endl;
		}
		else
		{
			
			if(bOutputDetail)
			{
				outFile<<"======================================================================"<<endl;
				outFile<<setw(16)<<"Lidar Const"<<setw(16)<<"Start Time"<<setw(16)<<"Average Mins"<<setw(16)<<"Height"<<endl;

				for(int i=0;i<highIndex-lowIndex+1;i++)
				{
					outFile<<setw(16)<<setprecision(6)<<list[i]<<setw(16)<<GetTimeString()<<setw(16)<<GetAverageMinString()<<setw(16)<<setprecision(4)<<IndexToHeight(i+lowIndex)<<endl;
				}
				outFile<<"---------------------------------------------------------------------"<<endl;
				outFile<<setw(16)<<"LC Average"<<setw(16)<<"Start Time"<<setw(16)<<"Average Mins"<<setw(16)<<"Height - Height"<<setw(8)<<"AOD"<<endl;
				outFile<<setw(16)<<setprecision(6)<<result<<setw(16)<<GetTimeString()<<setw(16)<<GetAverageMinString()<<setw(8)<<setprecision(4)<<tempLow<<setw(8)<<setprecision(4)<<tempHigh<<setw(8)<<setprecision(4)<<AOD<<endl;
			}
			else
			{
				outFile<<"---------------------------------------------------------------------"<<endl;
				outFile<<setw(16)<<"Lidar Const"<<setw(16)<<"Start Time"<<setw(16)<<"Average Mins"<<setw(16)<<"Height - Height"<<setw(8)<<"AOD"<<endl;
				outFile<<setw(16)<<setprecision(6)<<result<<setw(16)<<GetTimeString()<<setw(16)<<GetAverageMinString()<<setw(8)<<setprecision(4)<<tempLow<<setw(8)<<setprecision(4)<<tempHigh<<setw(8)<<setprecision(4)<<AOD<<endl;
			}

			outFile<<endl;
			outFile.close();
			delete []list;

		}
		return result;
	}
}
int OneLine::GetCloudBottomIndex(float newThreshold, float newThreshold2, float newThreshold3)
{
	if(bIsSigma)
	{
		cout<<"Error : Can't get Cloud Bottom from Sigma data... , returning -999"<<endl;
		return -1;
	}
	else
	{
		if(!bIsNormalized)
		{
			NormalizeData();
		}
	}
	
	int maxPossibleCloudHeightIndex=HeightToIndex(MAXPOSSIBLECLOUDHEIGHT); 

	for(int i=FIRSTDATAINDEX+1;i<binsNum && i<=maxPossibleCloudHeightIndex ;i++)
	{
		if( (data[i]>newThreshold)
			||(data[i]>newThreshold2 && (data[i]-data[i-1])>newThreshold3))
		{
			return (i-1);
		}
	}
	return -1;
}

int OneLine::GetMixlayerIndex(float newThreshold, float newThreshold2, float newThreshold3)
{
	if(bIsSigma)
	{
		cout<<"Error : Can't get Cloud Bottom from Sigma data... , returning -999"<<endl;
		return -1;
	}
	else
	{
		if(!bIsNormalized)
		{
			NormalizeData();
		}
	}
	int MLH_Index=0;

	int maxPossibleMLHeightIndex=HeightToIndex(MAXPOSSIBLEMIXLAYERHEIGHT); 
	for(int i=FIRSTDATAINDEX+1;i<binsNum && i<=maxPossibleMLHeightIndex;i++)
	{
		//Not finished yet...
	}

	return MLH_Index;

}

